Vortices in a thin film superconductor with a spherical geometry 



On 
On 

S3 
33 



M. J. W. Dodgson and M. A. Moore 

Theoretical Physics Group 
Department of Physics and Astronomy 
The University of Manchester, M13 9PL, UK 
(June 7, 1996) 

We report results from Monte Carlo simulations of a thin film superconductor in a spherical ge- 
ometry within the lowest Landau level approximation. We observe the absence of a phase transition 
to a low temperature vortex solid phase with these boundary conditions; the system remains in the 
vortex liquid phase for all accessible temperatures. The correlation lengths are measured for phase 
coherence and density modulation. Both lengths display identical temperature dependences, with 
an asymptotic scaling form consistent with a continuous zero temperature transition. This contrasts 
with the first order freezing transition which is seen in the alternative quasi-periodic boundary con- 
ditions. The high temperature perturbation theory and the ground states of the spherical system 
suggest that the thermodynamic limit of the spherical geometry is the same as that on the flat 
plane. We discuss the advantages and drawbacks of simulations with different geometries, and com- 
pare with current experimental conclusions. The effect of having a large scale inhomogeneity in the 
applied field is also considered. 
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I. INTRODUCTION 

The existence and nature of melting transitions in two- 
dimensional (2D) systems remains a topic of interest and 
uncertainty in both theory and experiments. In partic- 
ular, the melting of the vortex lattice in a thin film su- 
perconductor has received much attention. A continuous 
melting transition from an Abrikosov lattice to a vortex 
liquid via a phase of hexatic order, as in the Kosterlitz- 
Thoulcss-Halperin-Nelson- Young (KTHNY) theory, was 
first proposed over fifteen years ago.El A "transition" in 
the current- voltage response of thin film, superconduc- 
tors has been seen in several experiments. The vortex 
system becomes depinned above some transition temper- 
ature (the depinned vortices will dissipate energy as they 
move, leading to an Ohmic response) and the rise of the 
resistivity of the film is presumed to be caused by the 
melting of the vortex lattice. However, it seems that this 
transition temperature varies with the strength of the 
pinning centers in the thin films. For films of amorphous 
NbO the pinning is weak and no features can be detected 
in the current-voltage characteristics which may be at- 
tributed to a "transition" . This has been interpreted 
being due to the absence of the Abrikosov vortex lattice 

The superconducting state arises from macroscopic 
correlations in the phase of the Cooper pair wave func- 
tion, or equivalently the Ginzburg-Landau order parame- 
ter. This phase coherence is usually a consequence of off- 
diagonal long range order (ODLRO), i.e. non- vanishing 
(-0), where ip is the superconducting order parameter. It 
has been predicted that thermal excitations of the shear 
modes of the vortex lattice should destroy ODLRO in the 
crystalline state and the length scale itwer which phase 
coherence persists has been estimated. B'LJ 

O'Neill and Moore have performed Monte Carlo (MC) 
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simulations of thin film superconductors within the low- 
est Landau level (LLL) approximation using a spherk 
cal geometry to minimize finite size boundary effects.Q 
No freezing transition out of the vortex liquid phase was 
found in the temperature region investigated. The corre- 
lation function associated with the positions of the vor- 
tices was calculated and the associated length scale of 
these lattice-like correlations was found to diverge only 
in the zero temperature limit. It was noted that this 
happened on a similar length scale over which there was 
expected to be phase coherence from the arguments of 
Refs. |and0 

In the work of O'Neill and Moore, the superconducting 
film existed on the surface of a sphere, and this geometry 
is the one employed here. The magnetic field perpendic- 
ular to the surface is generated by .a, magnetic monopole 
placed at the center of the sphere.El The strength of the 
monopole must be quantized to have a well defined super- 
conducting order parameter, but the strength of the field 
at the surface may be continuously varied by changing 
the radius of the sphere. The LLL approximation, which 
assumes an infinite magnetic screening length, works well 
in the 2D case where the effective screening length is 
A e ff = 2X 2 /d because of the magnetic field outside the 
superconductor A e // becomes large as the film thick- 
ness d is reduced; A is the bulk magnetic penetration 
length. It is often the case that A e // is greater than the 
film dimensions for small d. 

More recently, strong evidence for a first order tran- 
sition at a finite temperature in an alternative finite 
size approximation with quasi-periodic boundary con- 
ditions (QPJ_an. an infinite plane has come from MC 
simulations £!j~hil This model has a periodic amplitude of 
the order parameter and a fixed change in the phase as 
a unit cell containing N vortices is crossed (the change 
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in phase depends on the choice of gauge). The quasi- 
periodicity in the phase places some constraints on the 
motion of the vortices, but it is not obvious that this 
will change the thermodynamic properties of the system. 
The observed transition temperature is within the range 
investigated on the sphere. This raises the question of 
whether or not there really is no transition on the sphere. 
If so, then what are the reasons that different boundary 
conditions give different results? 

If there is a freezing transition, at least on the plane, 
does this contradict the theoretical prediction of no 
ODLRO at finite temperatures? Apparently not, as the 
work of Sasik, Stroud and Tesanovictil (SST) has con- 
cluded that in the thermodynamic limit even the low 
temperature state below the first order transition has no 
ODLRO, despite the presence of "Bragg peaks" in the 
correlations of the order parameter density (the peaks 
will not be delta-functions but the algebraic singularities 
expected for a 2D crystalline systemt£l). These conclu- 
sions of SST are taken as support for the "charge-density 
wave" phase (with a modulation in the Coopec-pair den- 
sity, but no ODLRO) proposed by Tesanovic.113 

In this work we again fail to find a phase transition in 
the spherical geometry. We have investigated the correla- 
tions associated with both the phase and the amplitude of 
the superconducting order parameter. One advantage of 
the spherical geometry is that the liquid state may be in- 
vestigated over a larger range of temperatures. This has 
allowed us to compare the growing length scales of the 
phase and amplitude correlations over a large variation 
in temperature. 

Secti on [CJ] contains the formulation of our model. In 
section [III| it is shown how the high temperature pertur- 
bation theory on the spherical geometry gives the same 
free energy per vortex as the conventional theory in the 
thermodynamic limit. In Section IV we describe how the 



ground state of the vortex system on the sphere is af- 
fected by the necessity of containing twelve disclination 
defects within the triangular lattice. We have previously 
concluded that these defects with long range strains will 
give a finite energy cost per vortex relative to the flat 
plane ground state even in the thermodynamic limit .EJ 
However, the presence of dislocation defects were not con- 
sidered, and these can screen the strains of the disclina- 
tions to recover the thermodynamic limit of the triangu- 
lar lattice on a flat plane. This screening is described in 
the framework of elasticity theory, and the numerical ev- 
idence for dislo cations in the ground state is reexamined. 
In Section V A our calculations on the four-point density 
correlation function of the vortex liquid are shown. This 
is followed by o ur investigations on the phase coherence 



in Section VB 



The spherical geometry is in principle realizable using 
the end of a long thin solenoid to approximate the field 
from a monopole and placing it inside a spherical super- 
conducting film. However there would always be a prob- 
lem in the placement of the solenoid which would lead to 
a large scale inhomogeneity in the field strength at the 



spherical surface. Additional simulations are described 
in Section VI, where we add just such an inhomogeneous 
magnetic field. 

Finally, in the light of different conclusions from sim- 
ulations of the LLL system in different geometries, we 
disc uss s ome of the problems in these geometries in Sec- 
tion VII and comment on the experi ment al situation. 



Some conclusions are drawn in Section VIII 



II. DESCRIPTION OF THE MODEL 

We model a thin film superconductor by a spherical 
shell of radius R and thickness d <C R. An external 
magnetic field H(r) = rH(r), always directed perpen- 
dicular to the shell's surface is created by a magnetic 
monopole at the center of the sphere. We neglect all fluc- 
tuations in the magnetic induction so that it always has 
the mean value B — fiQH(R) on the superconducting sur- 
face. Dirac's quantization conditioru (assuming that the 
basic unit of charge is that of a Cooper pair) requires that 
the total flux through the surface is an integer multiple 
of flux quanta Nh/2e. Our choice of gauge that satisfies 
VxA = BisA = (A r , Ae,A^) = (0,0,B J Rtan6»/2). We 
also assume that the superconductor is described by a 
Ginzburg-Landau (GL) complex order parameter ip with 
a phenomenological free energy Hamiltonian given by 
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(2.1) 



where D 2 — D.D and D = — ih.'V — 2eA is the gauge 
invariant derivative operator. The LLL approximation 
consists in expanding this order parameter in terms of 
the eigenstates of D 2 , and then only keeping the degen- 
erate level of states with the lowest eigenvalue. If we only 
consider the two dimensions on the surface of the sphere 
then the subspace of-this LLL is spanned by N + 1 or- 
thonormal functions:t3 

1> m (9, <f) = hme"^ sm m (0/2) cos Ar - m (0/2), (2.2) 

with to = 0, ...A and h m = [(N + AnR 2 m\(N - to)!] 1/2 . 
We write the order parameter as ip(8, <fi) = 
QX>m "0m (#,</>) with Q = (tPokBT/fidB) 1 / 4 , and mea- 
sure lengths in units of the magnetic length l m = 
(^o/SttB) 1 / 2 which fixes R = {N/2) 1 ' 2 . [Note that 
the lattice spacing i n a triangular vortex lattice will 
be given by l = (4Tr/V3) 1/2 l m ^ 2.69l m .} The LLL 
Hamiltonian can then be written as 



H({u m }) = k B TctT 
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The quartic coupling term is 
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Wpqrs - 2N(2N + l) JpqJrs ' 
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2 - C p "Cf /C™ q and C\ = - is the 



where f pq 

binomial coefficient 

Our effective reduced temperature variable is 
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where a' = T o(da/dT)T=T c0 and T c q is the mean field 
transition temperature. H C 2(0) is the straight line ex- 
trapolation to zero temperature of the upper critical field, 
within mea n fie ld theory. For better-computational effi- 
ciency, Eq. (2.3) maybe rewritten asEa 



H ({u m }) = k B Ta.T 
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where U p = 2wNY jq=Q B x / 2 {2N - p + l,p + 
l)hqh p - q Q(p - q)Q(N + q - p)u q u p ^ q and B(x,y) = 
T(x)T(y) /T(x+y) is the Beta function. 8(<?) is the Heav- 
iside step function. 

An alternative description of the order parameter on 
the sphere in the LLL approximation is possible if we 
write (following Ref. 0) if>(0, 4>) = Q J2m=o v m^ m {0, 4>) = 
Q cos N (8/2) YZ=o v ™z m , where z = tan(0/2)e^. The 
Nth order polynomial may be written as a product of 
N + 1 terms to give: 



TV 



i/)(0, tf>) = Ccos N (6/2) ]J(z - 



(2.7) 



Therefore the order parameter is determined by the posi- 
tions of the N zeros (vortices) {zi} along with an overall 
complex amplitude C. It is shown in Appendix [X] that 
the GL free energy Hamiltonian depends only on the rel- 
ative positions of these vortices plus the overall complex 
amplitude, and not on the choice of coordinate system 
(6, (j>) that arose from our choice of gauge. This transla- 
tional invariance is assumed throughout this paper. Note 
that there are no constraints on the allowed positions {z{\ 
which is not the case with QP boundary conditions. 

The properties of these systems are de- 
termined by the partition function Z = 
I Um dumdu^ exp [-Ti({u m })/k B T] . We investigate 
these properties, in particular correlations in the "den- 
sity" \ip\ 2 and the phase arg(-0), with the Metropolis 
algorithm w-kuch uses an MC method to sample the 
phase space.El We have made measurements for runs up 
to 10 7 MC steps for system sizes as large as N = 400. 
We have mainly concentrated on the temperature range 
between ar = — 2 and —12. We have also used standard 
minimization routines to investigate the ground states of 
this model. 



III. HIGH TEMPERATURE PERTURBATION 
THEORY ON THE SPHERE 



A perturbation exp ans ion has been developed about 
the Gaussian limit of (2.1) at or — * +Q©ito approximate 
various properties of the LLL system.Ei In two dimen- 
sions the geometry used has always been a flat plane. 
The series may be represented by Feynman diagrams and 
many quantities in twOpand three dimensions have been 
evaluated to high order r3 The series can extend to nega- 
tive ut as long as the propagator is renormalized at least 
to the level of the Hartree-Fock approximation (HF). 

We consider the perturbation series for the free energy 
using the spherical geometry that leads to the Hamilto- 
nian (2.6). The usual treatment follows for the Feynman 
diagrams in a </> 4 theory, with directed lines representing 
the gaussian propagator go — (v p v*)o connecting vertices 
representing the interaction f3 pqrs = w pqrs 5 p + q , r + s , and 
internal lines are summed from to N. All diagrams with 
loops are removed by renormalizing the propagator in a 
self-consistent fashion. This is the HF approximation, 
and the renormalized propagator, g, must satisfy 



g = g (l-2g 2 S v ), 



where S° = X^fLo Ppqpq- O n the sphere we find: 



^sphere 



2N 



1 

N' 



(3.1) 



(3.2) 



so that in the large N limit on the sphere, the HF prop- 
agator approaches that on the flat plane. 

Exactly the same diagrams appear in the free energy 
expansion for the spherical geometry as those for the 
plane. The contribution of each diagram with n vertices 
to the coefficient of the nth order term of the expansion is 
a product of a combinatorial factor and the integral/sum 
over the "momenta" along each line of the interactions 
at each vertex. The only differences on the sphere are 
these sums. 

We now look at the lowest order Feynman graph with- 
out loops, shown in Fig. |l|a, which has two vertices con- 
nected by four lines. The sum that appears in this dia- 
gram is T 2 = J2 pqrs =o PpqrsPrspq- On the plane this is 
equal to T 



plane 



N/8, and we find on the sphere: 



T. 



sphere 



(N + l) 4 _ N J_ 
AN 2 {2N + 1) 8 ' + 2N 



■0(N- 2 )}, (3.3) 



which also gives the result on the plane with a 1/JV cor- 
rection. 

Generalizing our method to consider higher order 
terms, we calculate a particularly simple type of diagram 
to arbitrary order: the ring diagrams shown in Fig. [l]b. 
For such a ring diagram of order n, the relevant sum 
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We find 



(N + l) 



•111 
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sphere 2"A fn (2iV+l) 

= ^ a „ e [l + ^^+0(iV- 2 )], (3.4) 
where the result on a plane from Ref. E3^ is T 



plane 

2 1 ~ 2n N. So the contribution to the free energy to nth or- 
der approaches that on the plane as long as n is much less 
than the system size N. Although we have not verified 
this for more complex diagrams, we expect a similar re- 
sult. This is because the interactions at each vertex have 
only 1/N (and higher power) corrections on the sphere 
from the interactions on the plane, so with n vertices in a 
diagram, there should be a total correction of order n/N . 
Therefore as long as the perturbation theory is useful (ic. 
on the high temperature side of any phase transition) the 
thermodynamic limit on the sphere should be the same 
as that for the flat plane. Intuitively, it is also reason- 
able that the free energy per vortex in the liquid phase 
should be the same on the flat plane and on the sphere. 
In the liquid phase, correlation lengths are finite (see Sec- 
tion ^) , and so long as they remain small compared to 
the sphere radius, then the liquid on the sphere is "un- 
aware" of its curvature. 





(a) (b) 

FIG. 1. Simple Feynman diagrams in the Free Energy se- 
ries that we calculate on the spherical geometry: a) a second 
order diagram; b) part of the nth order ring diagram. 



IV. GROUND STATES ON A SPHERE 

The ground state of the LLL vortex system on an in- 
finite plane was long ago discovered to be the triangular 
lattice.O This is also the ground state on the finite sys- 
tems used in the QP simulations when commensurate 
boundary conditions are chosen. However, on a sphere 
the closest configuration the vortices can make to an ideal 
triangular lattice must contain twelve "disclinations" , i.e. 
twelve vortices that only have five nearest neighbors. 

We have previously studied the ground states by nu- 
merical minimization p£ the Hamiltonian in Eq. ( |2.6| ) for 
different system sizes JI3 The following conclusions were 



made: there exist some values of N (magic numbers) for 
which the ground state has a lower energy than nearby 
values, and the energies as N — > oo do not seem to con- 
verge to the infinite plane value. Extrapolation of our 
numerical results gives an extra energy on the sphere 
in this large N limit of 0.2(4)% of the ground state en- 
ergy per vortex on the infinite plane. The vortex con- 
figurations of these ground states make up a triangular 
network, but with twelve five-coordinated centers. The 
result that as N — > oo the ground states do not look the 
same as for an infinite plane may be explained by the 
presence of the disclinations (five-coordinated centers), 
which are defects with long range strains rB Therefore 
there will always be a finite fraction of the vortex system 
that is distorted from an ideal triangular lattice. This is 
supported by calculations within elasticity theory (which 
assumes only small distortions) that give a similar finite 
energy cost of the distortions on a sphere in the large N 
limit when an icpaahedral arrangement of the disclina- 
tions is assumed. lill 

This result of an extra energy per vortex in the ther- 
modynamic limit would be equally valid at finite temper- 
atures as long as there is a non-zero shear modulus (ie. in 
a solid phase). If it were true there could be a paradox: 
Imagine a 2D system with a continuous phase transition 
from a liquid to a crystalline solid phase. Above the 
transition the free energies per vortex have been shown 
in Section [II to be the same on the sphere and the plane 
and so this transition must remain at the same tempera- 
ture in each case. As at the transition the free energy is 
continuous, then the free energies in the solid state must 
also be equal. 

In fact, one can do a better job on the ground states on 
the sphere by having a number of dislocations (which are 
bound pairs of five-fold and seven-fold disclinations) in 
addition to the twelve free disclinations. This means we 
have to take away the restriction that every vortex is six- 
fold coordinated other than the twelve five-fold centers 
used in the elasticity calculation of Ref. [l7|. There is no 
topological problem with this as one can pull an opposite 
pair of dislocations out of a perfect lattice. That extra 
dislocation defects may lower the energy may be seen in 
two different ways. First, the interaction energy between 
a dislocation and a disclination may be negative, depend- 
ing on the direction of the dislocation and this interaction 
energy grows faster with system size than the self energy 
of the dislocation. Alternatively, consider a region on a 
sphere bounded by three great circles that intersect at 
neighboring disclinations. Even in the large N limit this 
region will be distorted from a triangle (all three of its 
angles will be 27r/5) and so a perfect triangular lattice in 
this region will be incommensurate with the boundaries, 
where a series of dislocations will occur. So we have a 
picture of a ground state in the large N limit made up of 
twenty regions of perfect lattice, but with sides contain- 
ing dislocations, and corners containing disclinations. 

There are still competing effects which discourage dis- 
locations, and which may explain why they were not seen 
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in Ref. [l7|. First there is the curvature of the sphere 
which favors the distortions of a disclination, but which 
decreases as 1/R 2 . Second, there is the positive self en- 
ergy of a dislocation, which is generally eptpected to in- 
crease as the logarithm of the system sizeES Third there 
is the interaction between two dislocations which is direc- 
tional, and approximately proportional to the logarithm 
of the distance between them. To estimate the density 
of dislocations needed to reduce the strain energy of the 
disclinations we use the framework of elasticity theory. 



A. Elasticity calculation of defects on a sphere 

The elasticity theory of twelve disclinations on a sphere 
is described in Rcf. [17j. The method is briefly explained 
here and then generalized to allow the inclusion of dis- 
locations. We consider small strains from an ideal lat- 
tice, Uij(r). These strains are related to a stress field 
by Hooke's law, cry = 2/iUij + \ukk$ij- For the LLL 
system the elastic constant A diverges (the system is 
incompressible) and the shear modulus /i is given by 
fi = 0.48fesTay/47r/3^Z^j. The stress tensor has zero 
divergence which allows the problem to be expressed in 
terms of the Airy stress function, x, defined by <7y = 
CikCjidkdiX- Five-fold disclinations are defined by the 
change in bond angle by 2n/6 when a path encircles one 
such defect. Dislocations, on the other hand, are defined 
by their Burgers vector, but they are also equivalent to 
a disclination dipole pair. The stress field corresponding 
to a particular configuration of these topological defects 
is given by the solutions to the biharmonic equation on 
a sphere: 



^ V 4 X = -^ + E 2 i S{r "'•>+£ ^h%S(v r d ) 

" s=l d 

(4.1) 



where Ko is the 2D Young's modulus which in the LLL 
system is given by Kq = 4/i. The free disclinations are 
labeled by s, and d labels the dislocations with Burgers 
vectors h d . In the large R limit where we neglect the en- 
ergy of bending the lattice over the spherical curvature, 
the energy of a given configuration is 



Fel = 



1 



2Kn 



d 2 r (V 2 x) 2 = 



1 



2Kn 



d 2 ra kk 2 , (4.2) 



where the integrals are over the surface of the sphere. 

We will need to use coordinate frames rotated so that 
the reference axis passes through the defect. For instance, 
if a disclination is at the position {0 S , <f> s ) in our original 
coordinate frame, a general point (9, <ft) is described in 
the rotated frame by the coordinates 9'(s) (the polar an- 
gle from the disclination) and <j>'(s) (the azimuthal angle 
in the rotated frame relative to the <fi = <fr s direction), 
where 



cos0'(s) = cos 9 cos 9 S + sin0sin0 s cos(0 — (f> s ) (4.3) 
sin#'(s) sin^'(s) = sin sin(</> — 4> s ) (4.4) 
sin0 (s) cos(f)'(s) = sm9 cos9 s cos(4> — 4> s ) — cos8sm0 s . (4.5) 

For the reference direction of the azimuthal angle <j>'{s) 
in the new frame also to be rotated, a further transfor- 
mation 4>'(s) — » 4>" {s) isjequired (as in the standard 
treatment of Euler anglescH) . We use this extra rotation 
for the dislocation defects, but for simplicity we do not 
define the transformat ion here. 

The solution of Eq. ( |4.l| ) in the absence of dislocations 
was found to be x = S s =i Xs where the sth disclination 
is at (0,, cb«) and 



V 2 x, 



' kk \ 



|^{ m i[l-cos0'( S )]+l}. 



(4.6) 



The elastic energy cost of the twelve disclinations is 
E\2 = R 2 J rfOE s ffcfe(^, 0)] 2 , which is proportional to 
the area of the system, and gives a finite energy cost per 
vortex. 

If we now consider additional dislocation defects, we 
may write X = J2 s Xs + J2d Xd where 



K 



Xd 



tfdjSiT - V d ). 



(4.7) 



By differentiating the solution for a disclination, we find: 
K l sm9'(d)sm4>"(d) 



O'fefcl 



47ri? l-cos0'(d) 



(4.8) 



where we take the size of the Burgers vector to be the 
lattice spacing l , and <j>"{d) is the azimuthal angle about 
the axis through the dislocation measured from the di- 
rection of the Burgers vector. 

Within elasticity theory, we can consider separately 
the contributions of the self energy of each defect, and 
the pairwise interaction energy between the defects. The 
self energy of a dislocation on a sphere is 



E, 



Kpl 2 
' 32tt 



1 + cos(a/R) 



21n-[l 

2 L 



cos(a/R)] 



E r . 



(4.9) 



where a is a cutoff of the order of the lattice spacing 
(within a of the dislocation, the linear elasticity theory 
must break down completely) and E core is the energy of 
the distortion inside this cutoff. As expected, this be- 
comes proportional to In (R/a) in the limit of large R, 
where: 



lim (E d ) 

H — ^oo 



M[ 2 Mjy.) 



2 In 2 - 11 + E r . 



(4.10) 
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In the presence of more than one dislocation the total 
elastic energy will be 

.12 N d 

E tot = — /dME^+E^) 2 ( 4 - n ) 

J 8=1 d=l 

12 N d 

= E 12 + N d E d + E E E sd(r s ,r d , b d ) 

3=1 d=\ 

N d N d 

+EE^( rd ' rd '' bd ' bd ')' ( 4 - 12 ) 

d=l d'>d 

which is a sum of the self energies plus pairwise inter- 
actions between every defect (the interactions between 
disclinations are included in E12). 

If we consider a disclination at (0 s ,<p s ) = (0,0) and 
a dislocation at (9 d ,<fi d ) with a Burgers vector pointing 
perpendicular to the 4> — 4>d hnc (such that the five-seven 
pair in the dislocation point towards the disclination), we 
find an interaction energy: 



E sd (r s ,r d ,b d ) = - 



KplpR 
12 



(1 - cos6> d ) 1 

— In - 1 - cos (9 d ) 

sin Vd 2 

(4.13) 



Allowing the dislocation to rotate by an angle 7 will give 
the above energy multiplied by a factor cos 7. 




e 



- 0.2- 

FIG. 2. The interaction energy _E^ os = (K l R/12)e sd be- 
tween a dislocation and twelve disclinations in an icosahedral 
configuration, where the dislocation is on a line between two 
neighboring disclinations, with a Burgers vector perpendicu- 
lar to the line. The polar angle 9o between two neighboring 
disclinations is defined in the text. The polar angle between 
the dislocation and the disclination is 6 as shown in the dia- 
gram, where the spots represent the positions of disclinations 
on the sphere. 



The interaction energy between two dislocations, one 
at (9d',4>d>) = (0,0), the other at (9 dl <j) d ), with oppo- 
site Burgers vectors perpendicular to the line joining the 
dislocations, is found to be 



E dd (r dl ,r dl b d \b d ) 



K l 2 

47T 



In i[l — cos9 d 



(l + COS0 d ) 



(4.14) 



Again as expected, the total energy of a dislocation pair, 
2E d + E ddl in the large R limit becomes proportional to 
the logarithm of the distance between the two opposite 
dislocations. 

Notice that the interaction energy between a disloca- 
tion and disclination increases more rapidly with R than 
the self energy and the dislocation-dislocation interaction 
energy. As E sd may be negative (depending on the orien- 
tation of the Burgers vector), there is an instability to the 
formation of dislocations, for large enough system size. 
If we consider the twelve disclinations to be fixed at the 
corners of an icosahedron, we find the energy landscape 
of a single dislocation to have minima along the lines 
that join the disclinations. The energy Ef d ° 8 (which is a 
sum of the interactions between the dislocation and the 
twelve disclinations) is shown in Fig. |as a function of 
position on one of these lines, with a Burgers vector that 
is perpendicular to the line. The position is denoted by 
the polar angle 9 from one disclination, with a neighbor- 
ing disclination located at 9 = 9o = 2tan~ 1 [(v / 5— l)/2]. 
By symmetry, the energy is zero exactly on a disclination 
and at the midpoint between two disclinations, and there 
is one negative stationary point. 

For an estimate of when dislocations should first start 
to appear in the ground states as R increases, we consider 
the energy to put two dislocations of opposite Burgers 
vectors on each line joining disclinations, at the positions 
that minimize E sd . Such a configuration has an elastic 
energy cost of 

Etot - £12 + K ll[-QMl{R/h) + 2.391n(i?/Z ) + 0.907] 



-60E r 



(4.15) 



where we assume a = Iq and have only taken into account 
the dislocation-dislocation interactions between pairs on 
the same line. If we neglect E core we find that E to t = E12 
when R = 6.05^o which corresponds to N ~ 530. How- 
ever, with an estimate of E core ~ KoIq/16tt, this critical 
value increases to N ~ 1000. The critical value of N 
is sensitive to the estimated values of the core energy, 
and the neglect of extra dislocation pair interactions. 
Also, we have neglected the local curvature, which for 
small systems will favor the presence of the disclination 
strains, and therefore disfavor dislocations that screen 
these strains. Our estimates are compared with numeri- 
cal results later in this section. 

Having established that dislocations may lower the to- 
tal elastic energy cost on the sphere for large enough 



G 



systems, it remains to show that large numbers of dislo- 
cations may screen out the strain field of the disclinations 
so as to cancel the R 2 dependence of the elastic energy. 




FIG. 3. (a) A disclination defect on a flat plane without 
dislocations, (b) An example of dislocation series screening a 
disclination on a flat plane. 



We consider series of dislocations along the lines be- 
tween the disclinations. The sum of the strains from 
each defect may be approximated at large distances by 
an integral. For such a series with equally spaced dislo- 
cations a distance cIq apart, the resulting strain fields are 
those from a positive and a negative disclination at each 
end of the line, but with a charge of size q s — 1/c (where 
q s = +1 for a five- fold disclination). This charge arises 
from the ratio of the distance between the five-seven pair 
in a dislocation to the distance between dislocations in 
the series. This is seen by writing the sum of the disloca- 
tion strain fields as an integral. A series of dislocations 



running from (0, 0) to some point (8, 0) will have a total 
strain field of 

^(0,0) ~il [ a$,(9,<p)d9' 
«o Jo 

= ~o£(M) - (4.16) 

because the strain field af, (9, <fi) of a dislocation at (9' , 0) 
with Burgers vector perpendicular to the = direc- 
tion is just the differential of the strain field a s e , (9, (j>) of 
a disclination: <jg, = lo[d<Tg, /d(R9')]. 

If we imagine such series with c = 5 between each 
disclination and the midpoints to neighboring disclina- 
tions, as is shown near one disclination in Fig. ^|, then the 
strain field of each disclination will be exactly screened 
away. However, there will be an effective disclination of 
charge q s = 2/5 at each midpoint [the midpoint makes 
a polar angle 9 m = 9^/2 with the disclination, so that 
tan# m = (y/E— l)/2]. What is needed is a variable spac- 
ing that starts from 5Zo at the disclination, but diverges 
at the midpoint. We may require that the spacing has 
the scaling form c(9') = 5 + (R/lo)g(9') with the con- 
dition that g(0) — (this form is chosen as any higher 
power of R would lead to a spacing that increases faster 
than the length of the series). 

In the large N limit an exact cancellation of the R 2 
terms in the energy is obtained. The strain field from 
one series starting at 9 = and lying along = is 
given by 

r e m 

a c (9,(t>)= 4,(9,4>)c- 1 (9')Rd9' 
Jo 

= K^.-Mc- 1 ^)]"" - | j'" aJ,(S,*)fc^f>< 

'-l< + WHnk^^-^ m * h (4 - 17 » 

where 1(9, 4>) = J®" 1 (dg~~ x /d9')ag,d9' , which remains 
of order one as the system size increases. In the 
large R limit we write the strain field of the zth 
series that is connected to the disclination j as 

- ^(1/5)^(9, ct>)+g' 1 (0 m )(lo/RH J . 9 JO 
(l /R)Ii,j(6,(l>) = -(l/5)<rj + (lo/R)fi,j where * runs 
from one to five. The elastic energy cost due to all the 
defects, E tot = E 12 + E sd + E dd is 

1 12 

3,3 — 1 

, / 12 12 5 ^ \ 

E * = ^ R2ejj ' " hR £ E / 
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which shows that all of the terms of order i? 2 cancel, leav- 
ing terms of order one which depend on the gradients of 
<7 -1 (#') and the value of g~ 1 (9 m )- Our approximation 
in simplifying the effective charge at the midpoint to be 
proportional to 1/R ignores a correction of order 1/R 2 
to this charge. As the energy of a disclination is propor- 
tional to the area of the system, the energy cost of this 
ignored extra charge is also of order one. Therefore, the 
production of these series of dislocations of variable spac- 
ing, along the lines between disclinations, is an effective 
way of screening the diverging strains of the disclinations. 
The ground state energy per vortex in the large N limit 
is the same as that for an infinite flat plane apart from a 
correction proportional to TV -1 / 2 due to the core energies 
of the dislocations. 



B. Numerical evidence for dislocations in ground 
states 

We now consider the numerical evidence for the above 
conclusions. It is sufficient to minimize the Hamiltonian 



of Eq. (2.6) by minimizing the Abrikosov ratio 



Ha 
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(4.19) 



The triangular lattice on the flat plane has the valueo 
0A,a — 1.1596. In Ref. [Tt], the ground states of system 
sizes up to N — 652 were found numerically, and it was 
suggested that there was a finite energy cost per vor- 
tex in the large N limit, consistent with the presence of 
twelve disclinations. The existence of magic numbers was 
observed and explained, although for large systems this 
only has a small effect on the total energy. 

We have extended this investigation to larger systems; 
this is motivated by the above elasticity results which 
suggest that above a critical size, dislocations will appear 
in the ground states and these will screen the strains of 
the disclinations to recover the large N limit ground state 
energy of the flat plane. In Fig. [| the sort of configura- 
tion is shown that we expect to find near each disclina- 
tion. This figure demonstrates just how large the spher- 
ical system must be to accomodate the dislocations and 
screen out the strains over a large region. 

In Fig. ^ our numerical results are shown for ground 
states of systems larger than N = 652 up to N = 2432. 



We also include results from Ref. [l?] for smaller system 
sizes. It is important to note that, as the systems get 
large, there is an increasing number of metastable states 
that the minimization routines may fall into. In each 
case we have taken the lowest energy after several hun- 
dred runs, but this is always just an upper bound on the 
true ground state energy. 
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FIG. 4. The Abrikosov ratio of numerically found ground 
states for different system sizes. Only magic number systems 
are shown. The straight line fit to N^ 1 ^ 2 is only on data 
for N < 1400. Above this value, dislocations are seen in the 
ground states and Pa falls under the fit. 

Plotting the Abrikosov ratio against N~ x / 2 the results 
for N < 1200 fall on a straight line that extrapolates to 
give a finite difference for the flat plane as N — > oo. How- 
ever, this trend changes above N = 1200, with energies 
falling below the straight line. The results are explained 
by the reasoning that once systems are large enough for 
dislocations to overcome their self energy and mutual re- 
pulsion, they may screen the strains of the disclinations 
and reduce the total ground state energy. In the nu- 
merical ground states with N > 1200 we observe disloca- 
tions near the disclinations with the correct orientation of 
Burgers vector (the distribution of the dislocations may 
be irregular due to the difficulties in locating the true 
global ground state). 

The elasticity calculation combined with this numeri- 
cal work is strong evidence that the ground state energy 
per vortex in the thermodynamic limit is the same on 
the sphere as on the flat plane. It should be mentioned 
however that for the system sizes we look at in MC simu- 
lations, the ground states have a higher energy per vortex 
by 0.5-1% of the ground state energy on an infinite plane. 
The jump in the energy at the transition observed in the 
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QP systems is of a similar order of magnitude (~ 1%) 
as the energy difference of the spherical ground states 
from the QP ground states. It must be a possibility that 
the raising of the ground state energy inhibits a freezing 
transition to the low temperature state, which could ex- 
plain the apparent absence of a phase transition in our 
simulations on the sphere. 



V. SIMULATIONS ON THE SPHERE 



temperatures we look at the function C™ = xT l~P 2 an< ^ 
then normalize this function by its infinite temperature 
limit Cj ]00 (also given in Appendix |b|). We average our 
measurements over m, so our investigations center on the 

function A(fc) = Y! m =-i C F /[( 2l + 1 ) C l,°°] with k = l / R - 
A(fc) is the same correlator that has been looked at in 
a plane geometry, by Hu and MacDonaldlij as well as by 
Yeo and Moore,E3but modified for a sphere. Examples of 
this function for different ay and N are shown in Fig. 0. 



We now go on to measure the correlations in the 2D 
vortex system using MC simulations. In the following, 
the system sizes we often use correspond to the magic 
numbers with particularly low energy ground states. 
This is because a phase transition would be most likely 
where the ground state is least frustrated. Also, where we 
start the simulations in the ground state, we have most 
confidence that our numerical ground states are the cor- 
rect ones in the magic number cases. However, we fail to 
observe commensurability effects in the MC simulations. 
The magic numbers are a property of the ground state, 
and do not affect the vortex liquid over the temperature 
range investigated. 



A. The Density Correlation Function 

O'Neill and Moore have investigated the correlation 
function associated with vortex positions .0 This method 
is unwieldy as it requires finding the roots of an Nth. 
order polynomial at every measurement. It is more con- 
venient to look at the four point density coccelation func- 
tion (briefly investigated by Lee and MooreEj-and related 
to the correlator studied on the QP systemal3) 



X(r- 



(|^(r)|>(r')| 2 ) 



<b/>(r)| 2 >(|^(r')n. (5.1) 



This is most revealing in reciprocal space. On a spher- 
ical surface this means that we need to decompose the 
function into spherical harmonics: 



= J dCtdn' X (r - v')Y{ 



Within the LLL subspace, this correlator can be written 
as thermal averages over the coefficients u p : 
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(5.3) 



with I™ l given in Appendix |B| and p equal to the mean 
density ~p = (|'0| 2 )- To make comparison with different 



We have concentrated our measurements on the form 
of the peak in A(fc) at the first reciprocal lattice vector. 
In our units for a triangular lattice with a unit cell of 
area <j>o/B this is given by |G| ~ 2.694. In the entire 
temperature range we look at (cut > — 12) we find that 
a Lorentzian always gives an excellent fit to this peak 
(see Fig. |jj). This is consistent with an exponential de- 
crease of the density correlations in real space. The asso- 
ciated correlation length £d is easily extracted from the 
width 8 of the peak at half its maximum, as £d oc 1/6. 
In Fig. ^ we plot 1/5 against ax for different system 
sizes, which shows how the correlation length grows as 
the temperature is lowered, and reveals when finite-size 
effects become important. Down to ax = — 10 the re- 
sults for the larger systems appear to be independent of 
system size. For ax < —8, finite size effects for N = 132 
make its values of <5 _1 fall beneath the trend of the larger 



systems. A linear fit at low temperatures is at least cor 
sistcnt 
that £d 



. 3I H 

sistcnt with the results from vortex-vortex correlations!!] 



\a T , as cM 



For small \ax\ corrections to this scaling form are ap- 
parent, so that the range of this fit is small. Neverthe- 
less, it is clear that there is no jump in the value of £d in 
the vicinity of the observed transition of QP simulations, 
or the predicted KTHNY transition. For the integral of 
A(fc) under a given peak to remain constant we would 
expect the general farm A(fe) ~ £,Df(£,D(k — G)) in the 
vicinity of the peak.Ea This implies that the peak height is 
proportional to the inverse width. Using the parameters 
from our Lorentzian fits, this relation is demonstrated in 
the inset of Fig. 0. 








FIG. 5. Density correlation function with N = 200 
for reduced temperatures = —2,-5 and —10 (the peak 
heights increase with lower temperatures). The inset 
shows the density correlation function at qt = —8 with 
N — 72, 200 and 400. The solid line is a Lorentzian fit when 
N = 400. 
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FIG. 6. Inverse width of the first peak in A(fc) with tem- 
perature. The inset shows the inverse width against the peak 
height at different temperatures. The triangles, squares and 
diamonds represent N — 132, 192 and 402 respectively. The 
straight line is a fit to the N = 402 points for \oit\ > 8. 



These results should be compared to those of Yeo and 
Moore where A(fc) was calculated within a high tempera- 
ture approximation that uses-a non-perturbative method 
to go to lower temperatures.E3 In this work similar zero 
temperature scaling is seen down to ax = —10, although 
the peak heights do not grow as fast. 

In contrast, Hu and MacDonald,E3 who calculate A(fc) 
with an MC simulation on a plane with periodic bound- 
ary conditions, see an order of magnitude jump in the 
peaks at a transition temperature of cut ~ ~ 9. Below 
this first order transition, A(k) shows "Bragg peaks" at 
the reciprocal lattice vectors. However, in the temper- 
ature range above this freezing, the correlation function 
appears to be quantitatively the same as on the sphere 
at the same temperature. This is as it should be as the 
liquid phase on the sphere has the same properties as 
the liquid on the flat plane as long as the system size is 
sufficiently large. 



B. Real Space Phase Correlations 

It is long range coherence in the phase of the or- 
der parameter (ODLRO) that gives rise to supercon- 
ductivity. Within the LLL approximation the phase 
can only be changed by the movement of the vortices, 
so one might expect a simple relation between vortex 
correlations, and phase coherence. Sasik, Stroud and 
TesanovicEfl have defined a phase correlation function as 
<r(r— r') = (ip* (r)ip(r')) /~p. From measurements of this in 
their MC simulation with QP boundary conditions they 
found that in the thermodynamic limit er(r) fell to zero 
over the order of a lattice constant even when in the solid 
phase characterized by long range density correlations. 

As our model is always in the liquid state it is clear that 
cr(r) will always be trivially short ranged simply due to 
the translations of the vortex system. However we are in- 
terested in the phase correlations between points moving 
in the same patch of correlated vortices. It is the length 
scale associated with this type of correlation that is physi- 
cally relevant. That is, when this length scale approaches 
the limiting scale of a given sample (eg. the distance be- 
tween pinning centers) it will become superconducting. 
To investigate this length scale we need to measure the 
correlations in the phase from a fixed vortex. We also 
need to fix the rotations about this vortex (these would 
also destroy phase coherence) . In practise we can do this 
by (a) fixing one vortex at a pole (8 — 0) by setting 
wo = 0, and (b) fixing the nearest six vortices by keeping 

1, 2, ... 6. {uj } are the coefficients for 



Hi 



(o) c . 
u) tor i 



the ground state with a fixed orientation. These condi- 
tions effectively fix the order parameter in a local region 
around the central fixed vortex. We can then measure 
the correlation between the phase at a point r' = (9' , cf>') 
and at r = (0, 0): 

a fix (r') = (rm(r')/\r(0)Mr>)\) flx . (5.4) 



10 



Note the difference between this and the correlation func- 
tion of SST. As well as the hxing of the vortices, our def- 
inition contains no part dependent on the amplitude of 
the order parameter. In terms of the MC variables, the 
phase at 8 = is given by ^(0, <f>)/\ip(0, <j>)\ = uie^/|«i| 
where <f> is the azimuthal angle of approach as 8 — ► 0. Fix- 
ing this angle <$> to be zero, the phase correlation function 
may be written as 



N 



^ u m h m sin 



72) cos N - m (8'/2)e im < t ' 



«i £m=i u mh m sin" 1 (0'/2) cos w -™(072)e 



im<f>' 



fix 



The reader may notice a problem in measuring the 
phase of a complex function on the surface of a sphere. 
If a path is followed that encircles one vortex, the phase 
of ip will change by 2n. However, an ambiguity exists in 
the definition of the inside or the outside of a closed loop 
in a closed surface and this could lead to an ill defined 
phase. In fact, it is our choice of gauge that places on 
one axis a string (or infinitely long and thin solenoid) ter- 
minating at the monopole that removes this ambiguity. 
This choice effectively makes the South pole (8 = tt) the 

"outside" of the system, and the North pole (8 = 0) the 
(5.& 

nter. This results in a singularity in the phase of the 
order parameter at 8 — tt. Encircling the South pole 



FIG. 7. The top diagrams show the phase correlation function \<jf ix (8, (j>)\ plotted in the complex z-plane where 
z = tan (0/2)e i0 for ar = —6,-9 and —11, with N = 72. The bottom row shows the mean density (|^i(r)| )fi X /p with 
the same fixed coefficients as for the phase correlations. Notice that one of the fixed vortices is a five-fold center. 



without apparently encircling any vortex results in a 
phase change of 2 N tt (this is the origin of the Dirac quan- 
tization of the charge on a monopole). This does not im- 
ply a singularity in the supercurrent, as the gauge part 
of the current will exactly cancel with this wave function 
contribution. 

Our measurements of <Jfi X {8, <f>) are shown in Fig. fus- 
ing the projection z = tan(#/2)e l< ^. For comparison, we 
also show the mean density (IVK 1 *)! ) fix with the same 
coefficients fixed. A growing range of phase coherence 
as the temperature is lowered is demonstrated, and coin- 
cides with the growth in density correlations. The length 
scale ^phaae is found by plotting cru x against R8 and fit- 
ting a gaussian over the region outside the fixed vortices, 
i.e. for R8 > a with a = 2.5. We find that ^ p hase increases 
proportionally to \cxt\ in the temperature range ax < — 7 
(see Fig. ||), with again corrections to scaling visible for 
ax ^ —7. The similarity between Figures |^ and [s] indi- 
cates that the length scales for phase and density coher- 
ence grow proportionally to each other as the tempera- 
ture is reduced. The length scales are unchanged within 
the measured errors for the different system sizes stud- 
ied. A relation between the growth of correlations in the 
density \ip\ 2 and in the phase of ip is clearly observed. 

The scaling of phase coherence £ p hase ~ |o!T^m was 
predicted by O'Neill and Moore from two distinct argu- 
ments. One of these involved considering the thermal ex- 
citation of small amplitude shear modes from the ground 
state lattice to define a phase coherence length within 
the crystalline phase. The second argument considered 
the free energy cost of phase fluctuations over a finite 
region of the vortex system. It is this second argument 
that is directly relevant to the liquid phase with a finite 
length scale of crystalline order, and we repeat and am- 
plify it here. The free energy cost of phase fluctuations 
over a region with crystalline correlations was shown to 
be F e ff = (d/2) f d 2 rc 6e lf ll (V 2 ± ^) 2 , where $ is the phase 



of the order parameter. From this result, the free energy 
cost of a phase change of tt over a region of linear extent 
£ will be F(£) ~ ceadlf n ^~ 2 . When such excitations pro- 
liferate, phase coherence will be destroyed. Hence, there 
will be phase coherence only over the length scale £ when 
the free energy F(£) is of order ksT. The cut dependence 
of the mean-field elastic shear modulus cqq is known in 
the large field limit, cqq ~ Qy^r/rfl^, with the simple 
result that £ ~ ja-r^m- 
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FIG. 8. The length scale associated with phase coherence 
extracted from aji x (B,(f)) with iV = 192 against temperature. 



We believe that there is another kind of phase fluc- 
tuation within the liquid that can be identified with a 
topological defect- a "braid" . Consider a ring of vortices 
of a given radius that move around the circumference of 
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the ring. (In three dimensions this vortex movement may 
take place along the field direction to form a topologically 
constrained defect- hence the name, braid. We have pre- 
viously calculated the energies of such defects. cB) Such 
a movement may be created by a simple combination of 
changes to the basis coefficients. If we make the transfor- 
mation u m — > u' m = u m e n for m = 0, 1, 2, . . . v and make 
no change to the remaining coefficients m = v + 1, . . . N 
then such a braid movement will occur in a ring of vor- 
tices that enclose approximately v vortices (see Fig. ||). 




FIG. 9. The path of vortices under the change 

Um — » u m = u m e 17 for m = 0,1,2,...!/ as 7 goes from to 
2n with v = 20. In this case the original state is the ground 
state for N — 72. The vortices are plotted in the complex 
z-plane where z = tan (6/2)e 11 *. 

If 7 = 2ir the vortices in the ring will each hop to a 
neighbors position, with a phase change in the order pa- 
rameter of 2ir inside the ring (i.e. no essential change). 
Excitations will involve values of 7 less than 2tt which 
increase the phase at the center of the ring by 7 while 
leaving the phase unchanged outside of the ring. The en- 
ergies of these braid excitations will define a length scale 
over which they proliferate, and thus destroy phase co- 
herence. The energy costs of braids in a crystal phase 
would increase with the size of the ring, and so phase co- 
herence would not be destroyed by this mechanism in the 
crystalline state. Nevertheless, we believe that the exci- 
tations of such defects are present on the length scale £ 
in the liquid phase, and are a plausible mechanism for 
the loss of phase coherence seen in Fig. [?]. Notice that 
in Fig. [7] it is when the density correlations become ring- 
like (as would happen due to braid excitations) that the 
phase coherence is lost. 



VI. EFFECT OF AN INHOMOGENEOUS 
MAGNETIC FIELD 

If an experiment using the spherical geometry was ever 
to be realized, there would always be a problem in plac- 



ing the end of a solenoid (the "monopole" ) exactly at the 
center of the sphere. This will cause the external field to 
be inhomogeneous over the surface of the sphere. Such 
inhomogeneities are generic in experiments with any ge- 
ometry, so it is worthwhile investigating these effects. 

We consider displacing the monopole from the center of 
the sphere by a distance d along the z-axis. It is straight- 
forward to show that for small e = dj R 7 this induces 
a change in the vector potential of 8 A = eB^RsmOcj). 
An inhomogeneity is introduced into the magnetic field 
at the surface which has a radial component of B r 
I ' 1 1 + 2ecos6*). Substituting this new field into the free 
energy gives an extra term to first order in e of 

» N 
m— 

where e = (dh 2 Q 2 /2m)e. 

Before we look at our MC results with this extra term, 
we consider its effect within the high-temperature pertur- 
bation theory. With the addition of ( |6.1| ) to the Hamil- 
tonian, the HF equation becomes 

Qra ,9m (1 ^ffm ^ ^ Op ^pmpin ) ? (^'^) 
P 

with g m = ax + e(2m — N) /N. The sum in this equation 
is no longer independent of to. However, in the large N 
limit, Ppmpm has the form, 

lim (f3 pmpm ) = — X -= exp[-(p - qf/A 2 }, (6.3) 

where A = [(2N -p-q)(p + q)/4:N} 1 / 2 . The width of the 
interaction is of order N 1 ' 2 over most of the sphere. If 
the variation in the renormalized propagator g p is smooth 
over the range of N 1 / 2 (which corresponds to changes in 
real space being smooth over the distance l m ) then the 
sum over p is dominated by contributions when p ~ to to 
S ive J2 P 9p(3pmpm — 9m/2- Therefore we have N + 1 in- 
dependent HF equations g m = g m (l — <?„)■ Each has the 
same solution, but with an m-dependent effective tem- 
perature: g m = (a m /2)[-l ± (1 + 4/a4) 1/2 ] with a m = 
aT+t(2m~N) /N. This will be true for any change in the 

Hamiltonian of the form SF = J2m=o e ( m )\ v m\ 2 ■ where 
the effective temperature becomes a m = aT+e(m). Such 
a change corresponds to changes in real space in one di- 
rection only over length scales longer than l m . 

This procedure can be generalized to an arbitrary dia- 
gram in the perturbative scheme. Any diagram is domi- 
nated by contributions where the indices in and out of a 
vertex are within A of each other. Therefore a given dia- 
gram is just the sum of the respective diagrams without 
the inhomogeneity but at the effective temperature a m . 
For example, the free energy may be expanded in a power 
series of the propagator, F pure (a T ) = NY^= Q a n g 2n . 
With the inhomogeneity this generalizes to 
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n—0 p—0 p—0 



(6.4) 



We also consider a third geometry with a circular geom- 
etry, where initial work has not found a phase transition. 
Finally, we look at the implications of recent experiments. 



This result may break down when variation in e(m) is on 
a scale less than A 1 / 2 . In the thermodynamic limit our 
conclusion will hold as long as the perturbation expan- 
sion remains useful. 
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FIG. 10. The m dependent propagator for two different 
strengths of inhomogeneity with N — 192 and e = —10. The 
dashed line represents the value of the propagator in the pure 
case at the effective temperature ar = ot m . Note that the 
effective temperature covers the range — 16 < a m < +4. 



Interesting effects may occur when there is a phase 
transition in the pure system. However, as our simula- 
tions fail to see such a transition, the above anal ysis is 
sufficient to explain our numerical results. In Fig. 10 we 
plot the m-dependent propagator |w TO | 2 from MC simu- 



lations with e = —10. and compare to the zero inhomo- 
geneity propagator with the m-dependent effective tem- 
perature. The results confirm that the system behaves 
with the local properties of the effective temperature. 



VII. DISCUSSION OF SIMULATIONS WITH 
DIFFERENT GEOMETRIES 



In this section we discuss some of the problems of dif- 
ferent boundary conditions when interpreting the results 
of numerical simulations of superconducting thin films 
within the LLL approximation. We first mention some 
of the shortcomings of the spherical geometry used in our 
simulations. We then go on to discuss the other main ge- 
ometry used- QP boundary conditions. In this case there 
exist problems not present in the spherical geometry due 
to the constraints on the phase of the order parameter. 



A. The spherical geometry 

An important difficulty on the sphere is that the crys- 
talline state free energy is raised due to the curvature 
and the necessary presence of twelve disclinations. A 
liquid with correlation length less than the system size 
will be unaffected by the curvature. The energy differ- 
ence should disappear in the thermodynamic limit (see 
Section [Tv| ). However, for the system sizes studied, this 
energy remains on the same order of magnitude of the 
jump in energy seen at the transition observed with QP 
boundary conditions. 

The presence of the topologically required disclinations 
on the sphere is clearly an important feature not present 
in real thin films. However, their effects will only be- 
come important as the correlation length of the system 
approaches the distance between disclinations. This is a 
finite size effect as experienced by any geometry, although 
the effective system size will be smaller for a given N than 
on a plane. For the largest system sizes studied, finite- 
size effects probably become important in the region of 
interest at ar ~ —10. 

The mobility of these disclinations, if large, would be 
a mechanism for destroying crystalline order in a finite- 
sized system. At low temperatures the disclinations 
should not be mobile, both because of their mutual re- 
pulsion (as can be seen from the elasticity calculation) 
and the topological barrier to their motion through the 
lattice. At higher temperatures where there is a finite 
density of dislocations (eg. in a hexatic phase) the dislo- 
cations provide a topological mechanism for the motion 
of the disclinations. In finite size systems, the strain 
field associated with the disclination may even create 
dislocations by pulling apart the dislocation pairs that 
are always present in a 2D crystal. From inspection of 
snapshots of the systems at finite temperatures we have 
observed that dislocations are attracted to the disclina- 
tions, and that the positions of the disclinations are not 
perfectly icosahedral as they are in the ground states. 

Our measurements so far have been insensitive to hex- 
atic bond-angle order. This order could never be globally 
defined on the sphere, and will only be a meaningful con- 
cept within the regions free of any disclinations, which are 
small in the system sizes we study. 

It must be conceded, therefore, that our numerical ev- 
idence from the sphere for the non-existence of a phase- 
transition is not yet conclusive. Larger system size stud- 
ies are clearly needed, but even increasing sizes by an 
order of magnitude may not resolve the problem. Our 
work on the ground states suggests that very large sys- 
tems with N ~ 10 4 may be needed to reach certain prop- 
erties of the thermodynamic limit such as the screening 
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of the disclinations by dislocations. 
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B. Quasi-periodic boundary conditions 
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In the system investigated by Refs. [To|-[l^ a constant 
magnetic field H perpendicular to a thin-film super- 
conductor is represented by the Landau gauge, A = 
—Byx. In this gauge the LLL states are tp p (x,y) = 
exp [ipkox — (y — pkol^ n )/2l^ n ], where fco is a "momen- 
tum" in the x direction that must be fixed by the bound- 
ary conditions. For a general state in the LLL sub- 
space, xp(x,y) — Q^2 p c p ip p (x,y), the system is peri- 
odic in the x direction with period L x = 2i: jk§. Quasi- 
periodicity in the y-direction may also be imposed by set- 
ting c m +N = c m . This condition leads to ip(x, y + L y ) = 
exp{iNk^x)%l){x,y) with L y — Nkol^. So each princi- 
ple region has an area L x L y = 2ttNI'^ 1 and contains N 
zeros of tp that correspond to vortices. To have a ge- 
ometry commensurate with a triangular vortex lattice, 
N x = L x /Iq and N y = L y /(^/3l /2) are fixed as inte- 
gers. (Note the difference between the magnetic length 
l m and the lattice spacing Iq both of which are defined in 
Section ||) 

The difference between this system and more conven- 
tional periodic boundary conditions comes from the con- 
straints on the phase of the order parameter as the prin- 
cipal region is crossed. That such a constraint exists is 
apparent from the fact that there are only 2N degrees 
of freedom available to describe the 2N coordinates of 
the vortex positions as well as an overall phase and am- 
plitude. It has been shown that the phase constraints 
lead ta the center of mass of the vortex positions to be 
fixed,Lj which accounts for the missing two parameters. 
A consequence of the fixed center of mass is that there are 
only N distinct ground states, each a triangular lattice 
separated by a minimum discrete distance. 

In MC simulations with QP boundary conditions, ev- 
idence is seen for a first order transition. The observed 
melting temperature drops as system size is increased but 
extrapolation gives an estimate for the thermodynamic 
transitior£3 at ar = —9.3 ±0.1. In the low tempera- 
ture state dominated by the ground state there is a finite 
energy barrier between different ground states. We now 
estimate what that energy barrier is. 



FIG. 11. Two possible ground states for N = 9 with QP 
boundary conditions. A path between the two ground states 
is shown that leaves the position of the center of mass un- 
changed. Note that only the principal region is shown. 



Figure |TT| illustrates the sort of movement of vor- 
tices between ground states that obeys the center of 
mass constraint. The example is for N = 9, and a 
short inspection will reveal that there are only 9 dis- 
tinct ground states that may be related by such moves. 
In general the different ground states have the form 
c p _ gi = exp [i(ir/2)(p/N x ) 2 - 2in(p + N xqi )q 2 /N] for p 
equal to an integer multiple of N x , and zero otherwise. 
Each ground state is labeled by ip qi tQ2 , where q\ must be 



an integer from 1 to N x , and the condition c p = c p+ n 
quantizes the allowed values of q 2 to integers from 1 to 
N y . These ,-gjround states are exactly equivalent to the 
Eilenbergercil states 0(r|ro), with ro = (xo,yo) where 
x o = ^lo/Ny and yo = qi{^/%/2)lo/N x . These are re- 
lated by 0(r|r o ) = exp [2m(2yo/y/3)x]<j)(r + r |0). 

Consider the change c p — > c' p = c p e %pl for any con- 
figuration of {c p }. For an infinite system this just cor- 
responds to a uniform translation of the vortices in the 
x direction by j/kQ. However, with QP boundary con- 
ditions this change will involve movement of vortices at 
the edges of the principal region in the opposite direc- 
tion, thus conserving the center of mass. (An analogous 
chan ge in the spherical geometry gives the braid in Sec- 
It is exactly this sort of movement that is 



tion 



VB 



a path between different ground states. For example, if 
we start from ^o,o & n d apply this transformation with 7 
changing from to 2irq2/N we arrive at the state ipo,q 2 - 
What is the energy barrier for such a path? The free 
energy of the QP system is given by 



W ({c m }) = k B Tc? T 



N-l 



sgn(a T ) I 

m=0 
I N-l 00 

2N~ ^ ^ 13 • 



*pqrs Cp Cq C r C s 



p—0 q,r.s—— 00 



, (7-1) 
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With/V, - 3 1 /4 2 -l/2 (5p+gir+se {-[( P - 9 ) 2 + (r- S ) 2 ]V3V4^}^ 

Changing the phase of the coefficients does not affect the 
quadratic energy. Also, the contributions to the quartic 
term f3 pqrs c p c q c*c* s will be zero when p, q, r, and s cor- 
respond to the same region (ie. they are all between 
and N — 1). As the width of the gaussian in the quartic 
interaction is of order N x , contributions will be negligi- 
ble unless these variables are within ~ N x of each other, 
and the edge of the principal region. Therefore the en- 
ergy barrier of this path will be of order N^l 2 - This N 
dependence has been observed numericallyEa 



The measurements by SST of this low temperature 
state are clearly not equilibrated as not all of the ground 
states are being sampled; there should always be a trans- 
lationally invariant density for any finite system. It is also 
a possibility that energy barriers between ground states 
affect the evidence for a first order phase transition. For 
instance, the free energy of the liquid will be overesti- 
mated if the simulation is not run for long enough (the 
standard time for these simulations has been ~ 10 6 MC 
steps) which may alter the melting temperature. This 
casts doubt on the extrapolation of the melting temper- 
ature to the thermodynamic limit. Even if a first order 
(or continuous) transition remains after full equilibration, 
the discrete ground state problem makes the QP system 
unsuitable for transport measurements. 



C. The circular geometry on a flat plane 



FIG. 12. The relative motion of vortices when moving 
between two neighboring ground states in the QP system. 

A more quantitative estimate may be obtained by con- 
sidering the movement of the vortices at the edge of the 
principal region relative to the vortices in the bulk of this 
region (see Fig |l2|). At any point on the path this repre- 
sents a defect from the perfect ground state, the energy 
cost of which may be calculated using the methods of 
Ref. ^9|. We find an energy barrier of 

E b = 0.33(L x /l )a 2 T , (7.2) 

which is proportional to N x as predicted. This is an up- 
per bound, as only the movement of a line of vortices 
has been allowed, without relaxation of the surround- 
ing lattice. However, as the interaction energy is short 
rangedjH3 we do not expect the N dependence to change. 
It is also a zero temperature result, but it should still 
hold for a low temperature state with long correlations. 
We emphasize that this energy barrier grows with system 
size, despite the fact that the positions of the discrete 
ground states get closer to each other. 

Additional evidence for the size of these barriers comes 
from the simulations themselves. In Ref. [l4| a plot is 
shown of (|V>(r)| 2 ) at a temperature corresponding to the 
extrapolated thermodynamic transition, where the mod- 
ulations in this quantity show that the system is dom- 
inated by just one of the N ground states during the 
length of the measurement, 10 6 MC steps. 



An alternative set of boundary conditions may be con- 
structed in the symmetric gaugeJ3 A = (xy — yx)B/2. 
The LLL states may be formed by the complete ba- 
sis with circular symmetry ip m — ( m exp[— |C| 2 /(4^m)], 
where C, — x + iy. To perform numerical simulations, 
one has to truncate the series of these states to a fi- 
nite number of terms, say, m — to N . This trun- 
cated system has N vortices, as tp = X)m=o c rai>m = 
exp [— |C| 2 /(4^)] rii=i(C — d): an< i the order parameter 
drops rapidly to zero outside of a radius {2N) 1 / 2 l m from 
the origin. 

Little work has been done on this geometry as finite 
size effects are strong. Numerical minimizations on small 
systems have shown that the ground state energy per vor- 
tex is lower than for the infinite plane,E3 and vortices near, 
the edge of the system are attracted to the "boundary"f!3 
(near where the order parameter is zero). The finite size 
effects are greater than with spherical or QP boundary 
conditions where there are no edges of the finite systems 
and each vortex resides in an equivalent local environ- 
ment. The effect of the boundary is mainly in a fraction 
of ~ N^ 1 / 2 of the bulk of the system, and should be 
unimportant in the thermodynamic limit. 

A similar geometry has been used in simulations of 
the 2D one component plasma (OGP), where the parti- 
cles are confined to a circular disk Ha A first order tran- 
sition was observed and this contrasts with the absence 
of a transitio£t,in simulations of the OCP on a spheri- 
cal geometry.E3 This has motivated us to perform sim- 
ulations on the circular geometry in search of a phase 
transition. We have performed the real space ph ase and 
density correlation measurements as in Section VB for 
systems up to N = 400 and find results consistent with 
those on the sphere. We have also measured the internal 
energy while heating from the ground state and cooling 
from high temperatures, for system sizes up to N — 1000. 
These energy measurements (see FigO) have revealed no 
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hysteresis on the scale observed in the QP system for tem- 
peratures greater than ax — — 12, and hence no evidence 
for a phase transition. 
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FIG. 13. The mean internal energy from MC measure- 
ments on the circular geometry with N — 1000. The cir- 
cles represent measurements while the temperature is being 
raised; the squares represent cooling. Each measurement is 
averaged over 30 000 MC steps. The energy is normalized by 
the ground state energy per vortex on the infinite plane. 



D. Experiments 

It seems that from a theoretical point of view there 
are three possible scenarios for the behavior of thin film 
superconductors in the presence of a magnetic field: a) 
the standard picture of KTHNY melting of a 2D crys- 
tal where a dislocation mediated continuous transition 
occurs between a low temperature state with quasi- 
crystalline order and a high temperature liquid or hexatic 
phase; b) a first order transition where the liquid free en- 
ergy becomes lower than the solid phase before KTHNY 
melting takes place (as is seen in QP simulations) ; c) the 
absence of a distinct low-temperature state, but where 
the correlation lengths of the liquid diverge as T — > 
(this is the prediction of high order perturbation theory 
from the high temperature side, and is consistent with 
the numerical simulations on the sphere). We now dis- 
cuss some experimental results and their implications for 
these possibilities. 

A melting transition of the flux lattice in In/In0 2 thin 
films was observed by Gammel et alH The zero field 
limit of the transition was consistent with the Kosterlitz- 
Thoulcss transition in the absence of external fields, 
though the magnetic field dependence of the melting line 
did not satisfy the predictions of the KTHNY theory. 
The measured transitions were mainly at low fields, so 
the relevance of our simulations here are doubtful. 

In the experiments of Berghuis et alB the current 
voltage response was measured of amorphous Nbi-^Ge^ 
films at different fields and temperatures. A drop to zero 
in the a.c. ohmic resistivity is seen, for example at fields 



H ~ 0.8H c2 at T = 0.68T C , though the d.c. response 
remains finite (ie. the critical current is zero.). This 
drop is attributed to a KTHNY transition consistent with 
the temperature and film thickness dependence of the 
melting field. A similar conclusion is made by Yazdani 
et alu from a.c. measurements on films of amorphous 
Mo 77 Ge 2 3. 

It should be emphasized that these a.c. measurements 
at a given frequency ui probe structural chapges of the 
vortex system over a restricted length scale.Q While this 
is useful to measure properties independent of the pres- 
ence of a large pinning length, one would expect even 
in the absence of a phase transition to see a large drop 
in the resistivity at fixed frequency when the relaxation 
times grow larger than 1/ui. In numerical simulations on 
the sphere,Ej the plastic time scale was observed to in- 
crease as T p i ~ exp (Aaj,) and we would expect a sharp 
drop when r p i ~ 1/uj. 

In contrast to the above experiments, measurements by 
Nikulov et ala on films of amorphous NbO^ with weak 
pinning have found no signatures of a resistive transition 
down to fields ~ 10 _3 £f C 2. This implies that the transi- 
tions observed previously are not universal properties of 
superconducting thin films. While the lack of a resistive 
drop may be attributed to the total absence of pinning, 
one would expect in the presence of a first order tran- 
sition, where there is a discontinuous change in (l^l 2 ), 
that this would be observed as a step in the flux flow 
resistivity. 

Thus the interpretation of these experiments remains 
controversial. While the KTHNY theory provides a use- 
ful framework to discuss flux lattice melting in two di- 
mensions, itsxffects have only been observed over short 
length scales .□ As far as we know there has been no ex- 
perimental indication of a first order transition in super- 
conducting thin films. In principle, the results of these 
experiments can still be explained simply by the growing 
length scales and relaxation times of a vortex liquid in 
the absence of a phase transition. 



VIII. CONCLUSIONS 

We have shown that the model with a spherical geom- 
etry remains in the vortex liquid phase, at least down to 
temperatures cut —12 for system sizes up to N < 402. 
This could be due to the raised energy of the ground 
states because of the necessary presence of twelve discli- 
nations in a triangular network over a sphere. However, 
we have made similar investigations on an alternative fi- 
nite size model: considering the LLL basis states on an 
infinite plane with a circular symmetric gauge and trun- 
cating the series to include only N vortices.13 This system 
has a lower ground state energy per vortex than the in- 
finite system and no disclinations, but we see the same 
length scales as in Section V B for temperatures down to 
ax = —12. The main evidence for a first order transition 
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in the QP model is the double peak structure observed in 
the probability distribution of the internal energy at the 
transition. Our investigations on the internal energy for 
both the sphere and the plane with circular coordinates 
have revealed no such signatures. 

The fact that the spherical system remains in the liq- 
uid phase down to low temperatures allows us to mea- 
sure the variation of the correlations in the liquid state 
with temperature over a large regime. Measurements in 
the liquid phase have the advantages that the system 
reaches equilibrium over a reasonable time scale at low 
temperatures. Our measurements of the density correla- 
tions are consistent with those in the liquid phase using 
QP boundary conditions Ji3 The extracted length scales 
demonstrate an asymptotic linear dependence on \ar\- 
The scaling analysis of O'Neill and Moore predicted that 
the phase coherence length should grow as t; p hase ~ \ a T\- 
Our results for the phase coherence are consistent with 
this scaling form, and also demonstrate that in the liquid 
the phase and density correlations grow with identical 
temperature dependence. 

Our results for the scaling of the phase correlations 
contradict the conclusions of SST, who claim that in the 
thermodynamic limit N — > oo the phase coherence has 
a range only of the order of the lattice spacing in their 
"solid" phase. However, this discrepancy is really due 
to different definitions of what is meant by phase co- 
herence. SST do not predict the growth of the length 
scale of phase correlations as the temperature goes to 
zero. Although their results are claimed, as evidence for 
the charge-density wave of Tesanovic,ll3 a finite length 
scale that grows with decreasing temperature or even a 
power law form of ODLRO would still be consistent with 
this exotic phase. It seems unlikely to us that the low- 
temperature phase below a first order transition would 
have a shorter range of phase coherence than the liquid 
phase. If we measure the phase correlations without fix- 
ing any vortices, we also find no phase coherence beyond 
a length of one lattice spacing, due to the translational 
freedom of the liquid phase. In the solid phase of the 
QP model the vortex system appears to be fixed near 
one ground state due to the finite energy barrier between 
the N ground states, and this allows phase coherence to 
be observed over a long range. As N increases, there is 
a decrease in the range of phase coherence. These mea- 
surements by SST demonstrate just how strong the finite- 
size effects are in the QP simulations below the observed 
transition. Without fixing the motion of any vortices, the 
phase coherence should always have a trivial length scale 
of the inter- vortex spacing. By restricting the motion of 
a small number of vortices, we have been able to measure 
a more physically relevant phase correlation length that 
is consistent with the expected growth of order as the 
temperature is reduced. An interesting question remains 
as to whether a similarly defined phase coherence in the 
solid phase with QP boundary conditions would result in 
true long range order, or another length scale that also 
exhibits zero temperature scaling. 



Despite the differences in results from spherical and 
QP boundary conditions, our results from high temper- 
ature perturbation theory and on the ground states sug- 
gest that the thermodynamic limit will be the same on 
the sphere as the infinite flat plane. This contradicts the 
earlier statement of Ref. [l7] where the effect of allowing 
additional dislocation defects into the ground state was 
not considered. Even so, the results on the ground states 
show how large system sizes must be for certain prop- 
erties to approach the thermodynamic limit. A definite 
drawback on the sphere is the small effective system size 
compared to the number of degrees of freedom. 

The spherical geometry is perhaps not only relevant for 
simulational purposes. It is possible that an experiment 
could be constructed with an amorphous type-II super- 
conductor sputtered over the surface of a sphere, and one 
end of a long and thin solenoid inserted into the center to 
approximate the field from a monopole. In contrast an 
experimental realization of the QP boundary conditions 
is not possible (these boundary conditions are primarily 
of interest in extrapolating to the thermodynamic limit 
of an infinite number of free vortices). We feel that the 
spherical geometry, with its property of translational in- 
variance, remains a useful system for numerical simula- 
tions, especially to measure properties of the unpinned 
vortex liquid. An investigation into dynamical proper- 
ties with a view to examining-.the length scales associ- 
ated with non-local resistivityo is planned, which would 
not be possible using QP boundary conditions due to 
the absence of full translational invariance of the vortex 
positions. 
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APPENDIX A: DEPENDENCE OF THE 
HAMILTONIAN ON VORTEX POSITIONS 

In this appendix we will show that the GL free energy 
Hamiltonian in Eq. ( |2.6| ) depends only on the positions 
of the zeros in the order parameter on the sphere. First, 
from Eq. ( p.7| ) we have 

N 

ip(6, cj)) = C J|[sin(6»/2)e^ - z, cos(0/2)] 
i=l 
JV 

= D JJ[oj sin((9/2)e # - h cos(0/2)] 
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APPENDIX B: DETAILS FOR DENSITY 
CORRELATIONS 



with D = CY\f =1 Zi/bi and bi = aiZi. If = when 
9 = 8i and cp = then we may choose: 



a, = cos(9 l /2)e~ i ^ /2 
h = sm(6 i /2)e i ^ 2 . 



(A2) 



Now, the GL-LLL free energy is a functional that de- 
pends only on the value of IV'I 2 = |-D| 2 | Yl?=i l/«| 2 a * eacn 
point on the sphere. If we can show that \fi(9, cj))\ 2 de- 
pends only on the relative position of the ith zero then 
the Hamiltonian will be invariant with respect to coordi- 
nate rotations of all the zeros: 



\fi< 



cos(6i/2) sin(0/2)e lW -^ )/2 



To hnd the reciprocal space density correlator from 
ther mal averages of the LLL basis state coefficients (see 
Eq. (5.3)) we need the integrals 



= -[l-cos6"(i)] 
= sin 2 [6»'(i)/2] 



where 6'(i) is the polar angle measured from the axis 
through the ith zero to the point (9,<f>), defined by an 
analogous expression to Eq. 4.3. This is independent 
of the coordinate frame and so we have shown that 
\ip(9, 4>)\ 2 depends only on the relative positions of the 
vortices plus the overall amplitude. 



T m — h h 

p,q,l ~ ""P I 



o Jo 



d9d(j} e l(q - p)<t > sin 8 sin p+9 (0/2) 



X cos 



2N-p-q 



(0/2)Yr(6,4>), (Bl) 



which may be written: 



Ip,q.l ~ 3p-m,q( 1) h p h p -. 

4:np\(N - p + to)! 
X (N + m+l)\m\ 3 2 



(2J + l)(l + m)! 



1/2 



- m)\ 

m — I, m + I + l,p + 1 
m + l,N + m + 2 



(B2) 



■sin(^/2) cos(6»/2) e - iW -^ )/2 



= — [1 — cos(9) cos(6*j) — 8111(6*) sin(6*i) cos(0 — <f>i)\ 



where 3F2 [a, b, c; d, e] is a generalized hypergeometric se- 
ries. To find the infinite temperature limit of C™ we 
consider the case to = and calculate the averages 
(v p v*v r v*) and (v p v*) exactly in the limit oct — > +00 
(this involves only Gaussian integrals). Our result is 



(A3) 



C, 



{21 + 1) 



N 



(TV + 1)2 ^ 

y > p=0 

(AH) 2 



-/,/ + l,p+l 
l,N + 2 



(N-l)\(N + l + l)\ 



(B3) 



(We are indebted to T. Blum for this final result.) This 
is consistent with the fact that in the limit N — > 00, the 
infinite temperature correlation function should be iden- 
tical wi th t hat for an infinite planeji 2 ] anc j the expression 
in Eq. (p|) should tend to exp {-l 2 /N). 
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